Spatially offset optical coherence tomography: Leveraging multiple scattering for high-contrast imaging at depth in turbid media

The penetration depth of optical coherence tomography (OCT) reaches well beyond conventional microscopy; however, signal reduction with depth leads to rapid degradation of the signal below the noise level. The pursuit of imaging at depth has been largely approached by extinguishing multiple scattering. However, in OCT, multiple scattering substantially contributes to image formation at depth. Here, we investigate the role of multiple scattering in OCT image contrast and postulate that, in OCT, multiple scattering can enhance image contrast at depth. We introduce an original geometry that completely decouples the incident and collection fields by introducing a spatial offset between them, leading to preferential collection of multiply scattered light. A wave optics–based theoretical framework supports our experimentally demonstrated improvement in contrast. The effective signal attenuation can be reduced by more than 24 decibels. Notably, a ninefold enhancement in image contrast at depth is observed in scattering biological samples. This geometry enables a powerful capacity to dynamically tune for contrast at depth.

Here, is a factor converting power to current, 1 and 2 are two-dimensional vectors transverse to the optical axis, and the subscripts R and S refer to the reference and sample arms, respectively. | ( )| is the modulus of the normalized temporal coherence function of the source, which has a value of 1 when the sample and reference arm pathlengths are matched.
The field backscattered from the discontinuity in the sample arm can be expressed using the extended Huygens-Fresnel Green's function, ( , ), which describes the propagation of a wave through an arbitrary material with random inhomogeneities using ABCD matrix transfer formalism (19). The field in the mixing plane can be written as: ( ; ) = � ( ; ) ( , ) = � ( ; ) 0 ( , ) ( , ) where is the backscattered field in the discontinuity plane and r and p are the transverse vectors in the mixing and discontinuity planes, respectively. The Green's function can be further decomposed into the propagation in a system without inhomogeneities, 0 , and a random phase distorting factor (caused by the inhomogeneities). We assume that the statistical properties of the incident and backscattered fields are independent and that the statistical properties of the bulk tissue and the tissue discontinuity are also independent. Substituting Eq. S3 into Eq. S2 and omitting the implicit dependance for simplicity, we can write the sample arm mutual coherence function as: The first term in equation S4 leads to the mean backscattered intensity in the discontinuity plane, 〈 ( )〉. Assuming diffuse backscattering from the discontinuity (10): were, is the Dirac delta function, is the incident power in the sample arm, and is the Fresnel reflectivity of the discontinuity. Absorption is included through the introduction of the complex refractive index, = − , where the imaginary part of the refractive index depends on the absorption coefficient, , and the wavevector, : = 2 , where = 2 0 and 0 is the mean free-space wavelength. The first term in the brackets can be interpreted to represent the ballistically scattered light in the absence of tissue inhomogeneities while the second term represents the multiply scattered light. Here, we introduce the lateral offset in the beam, , in order to obtain the spatially offset backscattered intensity in the discontinuity plane: The quantities 2 and 2 are the 1 � intensity radii in the absence of scattering or absorption and in the presence of both scattering and absorption, respectively. These are given by (10,21,44): In these equations, and are the ABCD matrix elements and 0 is the 1 � radius in the lens plane. For our geometry, = 1 and = + (in the presence of absorption, B is complex). 0 is the lateral coherence length (19,46), and is given by: (1 + ( ) ). Also, (S10) The dependance of on z in Eq. S10 is indicated to include OCT systems that employ dynamic focusing, where = − .
The second term in Eq. S4 is the Huygens-Fresnel Green's function describing propagation from the discontinuity plane to the mixing plane in the absence of tissue inhomogeneities. The Green's function will differ from that presented in Thrane et al. (10) due to the incorporation of absorption. With this change, the Green's function becomes (47): where, 0 is the total optical path length from the discontinuity plane to the mixing plane, and , , and are elements of the ABCD transfer matrix. For our geometry, = = 1 and = + .
The third term in Eq. S4 represents the impact of tissue inhomogeneities on wavefronts propagating through the tissue, and can be expressed as the mutual coherence function, of a point source located at the discontinuity plane and observed in the lens plane (10): (S14) Finally, the electric field in the reference arm (ignoring complex terms that cancel when the coherence function is calculated) can be written as: where is the incident power in the reference arm. The presence of the offset, , in the reference arm depends on the geometry of the OCT system. If the offset in the sample arm can be adjusted without impacting the reference arm, such as the 800 nm system demonstrated in this study, the offset in the reference arm can be set to 0. We include it here for generality since it will often be relevant for Michelson-interferometer-based SO-OCT.
Substituting the relevant equations into Eq. S1 and performing the integrations over and , we now obtain an expression for the mean square heterodyne signal current: where is the mean square heterodyne current in the absence of scattering and absorption, and ( , ) is the heterodyne efficiency factor. ( , ) represents the tissue attenuation due to absorption and scattering, and can be expressed as: (S17)

Extension of EHF modelling framework to encompass angular offsets
We acknowledge that the implementation where the illumination and collection paths are offset and parallel represents one case of a broader framework encompassing the spatial and angular offset between the illumination and collection paths. Thus, we can extend our model to also incorporate an angular offset and connect our framework to the dual axis scheme (15). For an overview of the sample arm geometry incorporating angular offset with and without a lateral offset, we refer to Fig. 1c and 1d.
Since no integrations are performed along the axial depth, , Eq. S17 can easily be generalized to encompass an angular offset between the illumination and collection paths by incorporating a depth-dependent variable offset. The depth dependent offset can be expressed as ( ) = tan + 0 , where 0 is the lateral offset in the discontinuity plane. In the limit of the paraxial approximation, ( ) = + 0 . Thus, the final expression for the heterodyne efficiency factor incorporating both lateral and angular offset can be expressed as: (S18) Building on the work presented in (40), we note a key difference between spatially-offset OCT and dual-axis OCT. When α ≠ 0, some ballistic light scattered at a particular angular range is collected by the detector. Indeed, this configuration selectively collects ballistically scattered light from the elliptical region where the illumination and collection paths overlap. Conversely, if α = 0 and is greater than the lateral resolution, all ballistic light is rejected. The angular relation between the illumination and collection paths in dual-axis OCT is designed to selectively collect light scattered ballistically at a particular angle to the illumination beam and is based on the principle that multiply scattered light, while collected by the coherence gate, contributes only an unwanted background signal that results on the degradation of OCT image contrast. The impact of angular offset on the heterodyne efficiency factor (and thereby the OCT signal intensity) is shown in Fig. S1. Conversely, recent work modelling the OCT signal using the extended Huygens-Fresnel model has shown that multiply scattered light does indeed contain information about sample structure (10). As such, the spatial offset allows selective collection of the multiply scattered light, which enables deeper imaging since the multiply scattered light has typically travelled deeper into the sample. The relative dynamic range of images with lateral and angular offset can be seen by comparing the slopes of the lines plotted in Fig. 2d and S1a; a shallower slope indicates that a smaller detector dynamic range is required to visualize the signal. Note that in Fig. S1a, the normalized contrast would be equivalent to the heterodyne efficiency factor since the peak contrast is always at the sample surface. The optimal configuration including the lateral and angular offset must be determined based on the specific geometry and scattering properties of the sample. Since the offset can easily be tuned from s= 0 to various offsets with s>0, images with and without the ballistic component can be acquired using the same OCT system. The magnitude of the spatial offset and the angular offset represent two factors that can be independently tuned based on the sample scattering properties in order to optimize the contrast at different regions by controlling the ratio of ballistically and multiply scattered light that are collected from different regions of the sample.

Enhanced contrast in microbead phantoms: microparticles in PDMS gel
In order to further explore the contrast dependance on the angle-dependent scattering properties of the sample, we performed a second experiment to demonstrate how the relative contrast between different sized structures changes with increased offset. Using the OCT system with central wavelength of 1295 nm described in Fig. 2b we employed a phantom comprised of 0.1% w/w Ti02 microparticles with mean diameter < 5 µm uniformly dispersed in a PDMS gel. A ~1mm slab of this phantom was placed on top of a flat piece of transparent epoxy. It is clear from Figs. S2a-d that the bulk scattering contrast of the Ti02 microparticles, which fall within the predominantly forward and backscattering regime, is greatly reduced as the offset is increased, while the contrast of the bottom layer of the phantom relative to the bulk scattering is increased. To quantify this improvement, we calculated the contrast-to-noise ratio (CNR) (30) for the depth range between the blue lines (signal) and the depth range between the red lines (background "noise") for each B-scan. The full image width over the indicated depth range was used to calculate the CNR. The red region was selected to show a region of the image where the intensity of the smaller particles is reduced relative to the larger surface, indicated by the blue region. Here CNR is defined as | − | � 2 + 2 ⁄ , where , , , and denote the mean and standard deviation of the pixel intensities in the signal region and background region, respectively. We measured values of 1.01 dB, 1.29 dB, 1.43 dB, and 1.47 dB for offsets of s = 0, 30, 40, 50 µm, respectively, which demonstrates that the CNR at depth increases as the offset increases.

OCT imaging in soft tissue: krill eyes
We take OCT images from an adult krill eye to show the capability of our SO-OCT system to reveal features at depth (39). Figs. 6 c and d show the non-offset OCT and SO-OCT B-scans (s = 150 µm), respectively, on the eye of an adult krill. The top surface of the eye is well-defined in both B scans. However, the bottom surface of the eye is poorly defined because of rapid signal attenuation. The signals from the bottom surface can be enhanced when an offset, in this case of s = 150 µm, is applied. As shown in our modelling (Fig. 3), the selected offset specifies a depth in the image which has optimal contrast enhancement. Thus, the large offset employed here was selected to demonstrate an improved contrast deeper in the sample, at the depth corresponding to the bottom surface of the eye. A higher CNR value of 5.93 dB was obtained in SO-OCT compared to 3.35 dB in conventional OCT, which corresponds to an enhancement factor of 1.77 in the obtained image quality.

Detailed descriptions of experimental setups for SO-OCT
The first SO-OCT system ( Fig. 2A) is entirely home built and is based on a Mach-Zehnder interferometer. In this configuration, we use a superluminescent diode (SLD) laser (S850, Superlum, Cork, Ireland) with central wavelength of 800 nm and a bandwidth of 14 nm as our light source, which is split into two paths by a beamsplitter B1, 90%, into the sample arm and 10% into the reference arm. A low NA microscope objective (LSM03-BB, Thorlabs, Newton, NJ, USA) is used for focusing the light up to 1 mm into the sample. The focus plane is imaged on a CCD camera (CCD1), which is on the direction of zeroth diffraction order of the grating (GR1). The first diffraction order is focused onto a line CCD camera (CCD2), which provides the spectra information. The reference arm propagates through the same optical path except the dispersion compensation (DC) before it interferes with the sample arm. As indicated by the double arrow in Fig. 2A, the tuning of the spatial offset s is realized by simply shifting the lens L2 using a translation stage. The linear relationship between the shift of L2 and the spatial offset s is calibrated carefully in advance by measuring the laser spot on CCD1. In this way, it is straightforward to acquire OCT images with or without spatial offset in order to compare the performance. The exposure time of the camera was adjusted to ensure comparable dynamic range in all images.
The second experimental setup is an adaptation of a commercially available frequency domain OCT system (TELESTO-II, Thorlabs Inc., Newton, NJ, USA), with a custom-built add-on enabling spatially offset detection (Fig. 2B). This system is based on a Michelson interferometer and uses an SLD with a central wavelength of 1295 nm and a bandwidth of 217 nm. A circulator (CIR) was added in between the light source and the interferometer. A 50/50 beamsplitter (BS) splits the light into the sample and reference arms. A low NA microscope objective (LSM03, Thorlabs Inc., Newton, NJ, USA) with nominal resolution of 13 μm in air was used to focus light onto the sample. A pinhole on a translation stage was introduced in the collection path in order to tune the offset, as indicated by the arrows in Fig. 2B. One difference compared with the 800 nm setup is that, by placing the pinhole along the collection path, an offset is also introduced to the light collected from the reference arm. Thus, the power in the reference arm must be increased to compensate as the offset is increased. This could impact the sensitivity of the measurement especially for large offsets when the maximum reference arm power has been reached. All images were acquired in high sensitivity mode with an A-scan rate of 5.5 kHz. for various angular offsets. The circles plotted on each line indicate the depth with highest contrast. The following system parameters representing the system described in Fig. S4 were used for these calculations: = ; = . ; = ; = . ; = − ; = − ; = . . All angles are in units of radians. Note that the calculations here all consider an OCT system with dynamic focusing.